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1 ABSTRACT 



Glauber dynamics, applied to the one-dimensional Ising model, provides a tractable model for 
the study of non-equilibrium, many-body processes driven by a heat bath. We explain Glauber's 
dynamical Ising model in the context of a system comprising three "spins" in a 1-dimensional Ising 
model. We first present the model without a magnetic field and then, following Glauber, add an 
oscillating magnetic field. We then show that Glauber's demonstration of the fluctuation-dissipation 
theorem carries over to the 3-spin case. A web- link to a simulation code and an interactive video 
display of the "spins" flipping from a non-equilibrium to an equilibrium conflguration (zero magnetic 
fleld) is included. 



The study of thermodynamics traditionally deals with equilibrium properties of matter. However, as 
remarked by Le Bellac, et al.{JV\,p. 335) "this is rather restrictive since non-equilibrium phenomena, 
such as heat conduction or diffusion, are of great interest and cannot be ignored." On the other 
hand, our understanding of the properties of out-of-equilibrium matter is still something of a work 
in progress, as one might conclude from Kubo's brief summary of modern developments in the 
theory of irreversible processes. (j2J, p. 374). 

The Ising model p| considers arrays of coupled two-component spins, and is often discussed 
in connection with ferromagnetism in modern textbooks. The model lived most of its life as 
an equilibrium model, until, as discussed by Racz [4], Glauber [5] initiated an industry of kinetic 
spin models (currently, "Glauber Ising" gets about 1.4 million hits on Google) by providing for 
single spin-flips induced by a heat bath followed by a relaxation to equilibrium. A solution for 
an arbitrary number of spins was subsequently given by Felderhof [6| . Kawasaki quickly followed 
Glauber's paper with a spin-exchange model that conserved total spins. Garrido, et al. [8j, using 
Glauber dynamics with competing heat baths, demonstrated the existence of non-equilibrium phase 
transitions in 1-dimension with short-range forces. Lage pj extended the use of Glauber dynamics 



2 INTRODUCTION 



2 



to the Potts model, taking care to point out that the extension is not unique, that is, that there 
are multiple possible dynamics that could be added to the original Potts model. 

Many authors have recognized that the Ising model can also serve as a lattice-gas model, the 
two "spin" orientations being considered as lattice occupation numbers, and 1. Chomaz et a/.[10j 
invoke such an interpretation to discuss the formation and evolution of compact stars. The Ising 
model has also found application to the discussion of the thermodynamic properties of Schwarzschild 
black holes [T^ A recent paper by Mazilu and Williams "introduces nonequilibrium statistical 
mechanics" by solving exactly a two-temperature (1-d) Ising model with four spins. 

The Ising Hamiltonian for a lattice with spins is [13] 

N N 

H = -J^SiSi+i - H^Si (1) 

i=l i=l 

with index + 1 = 1. H is the strength of an applied external magnetic field, which may be time 
dependent. The coupling strength J is positive for an attractive interaction that would simulate 
ferromagnetism if all the spins had the same value, thereby indicating alignment. 

Statistical mechanics typically deals with very large systems, see e.g. [13]. Glauber's [5] paper 
accordingly considers the infinite 1-d Ising chain (as well as the single spin case). Chemists and 
biologists, however, are often interested in very small systems, such as molecules with just a few 
binding sites. Ter Bush and Thompson |14] worked out the Glauber dynamics of 2-spin and 4-spin 
Ising chains in connection with a model for the functioning of hemoglobin. 



3 GLAUBER-ISING WITH 3 SPINS 

The 3-spin 1-d dynamic Ising model displays many of the features of larger finite Ising systems 
and is easily handled by elementary methods in the case of zero magnetic field. It is assumed 
that each individual spin, Sj fiips randomly between the values +1,-1 with a probability per unit 
time (measured in arbitrary time units) that depends upon the orientations of the adjacent spins 
according to 

Wi{si) = [l- -7Si(si_i + Si+i)] (2) 

The constant 7 is related to the spin-coupling J and the temperature of a presumed heat bath, 
that drives the random orientation transitions of the individual spins, by [5j 

7 = coth(^) (3) 

where T is the temperature of the heat bath. Appendix A includes a table of the values of the ifj's. 

Designating the probability of a given spin configuration by p(si, S2, S3; i), where each Sj can 
take either of the values +, — , the time dependence of the probabilities is given by the master 
equation 

^^'^^^'^r^^'^'^ = - Y^Wi{si)p{si,S2,S3;t) + Y^Wj{-Sj)p{sj-i,-Sj,Sj+i;t) (4) 
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with the convention that subscripts and 4 are respectively the same as 3 and 1. As pointed 
out in connection with Eq. 17 of |^5j, the couphng strength 7 is hmited to the range [+1, —1], the 
extreme values corresponding to zero temperature of a presumed heat bath, and the negative value 
corresponding to a repulsive interaction between adjacent spins. The coupling vanishes in the limit 
of infinite temperature of the heat bath. 

Now let V-^{V-) represent the probability for the spins all having the value 1), and X^{X-) 
the sum of the probabilities that exactly two spins have the values +1(— 1). The individual prob- 
abilities for the latter configuration are denoted by x±i with the "i" denoting the position of the 
odd spin. Also, define V = V+ + V^, X = X+ + X_ , U = V+ -V_, Y = X+ - X_, Xi = + x_i, 
and Ui = x+i — x-i. Clearly, in the case of three spins 

V + X = l (5) 



4 ZERO MAGNETIC FIELD 



Define a vector 



V 
X 

u 

Y 



(6) 



Then, In Glauber's "alternative" approach, the master equation may be written as (the dot denotes 
time differentiation) 

4-^' = if = JVW (7) 



dt 



where the matrix M is 

-3(1-7) (1+7) 
3(1-7) -(1 + 7) 










-3(1 
3(1- 



-7) 
7) 






(1 + 7) 
-(5+7) _ 

?he eigenvalues of the matrix M are 0, —2(1 — 7), ~2(2 — 7), —6. The zero eigenvalue corresponds 
to equilibrium. There is a "parity" symmetry in the equations of motion, corresponding to the 
simultaneous flipping of all three spins, so that both U and Y must be equal to zero at equilibrium, 
and the equilibrium solution is ^ = 2(2'-^) ' ^ would be expected from the static Ising model. It is 



evident that the six configurations with one odd spin must be equally populated. 
The solutions are 



X{t) 
U{t) 
Y{t) 



X(0)e" 



-2(2-7)t , ^0--^) r. _ p-2(2-7)tl 
2(2-7)^ J 



^{[7(0)[(1 - 7)6-*^* + 3(1 + 7)e-2(i-T)* - [(1 + 7)y(0)[e-6* - e-2(i~T)*]} 



2(2 + 7) 
1 

2(2 + 7) 



{y(0)[3(l + 7)e-6* + (1 - 7)e-2(i-7)t] _ 3(1 _ 7)C/(0)[e 



-6i _ g-2(l-7)t 



]}(8) 
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The first of these equations provides for relaxation to the classical Ising equilibrium solution, as 
shown in the accompanying Figure, where, in the initial configuration the three spins are aligned 
in the positive direction (-'^(O) = 0, U{0) = +1, so that 1^(0) = 0). The scale of the t— variable is 
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4.1 Average Values 

Glauber also suggests that in many situations the relevant observable quantities might be average 
values and correlations of the individual spins, rather than the detailed configurations. Accordingly, 
he defines the single spin average value to be, for the k'th spin 

Qkit) =< Sk >=^Skp{si,S2,S3;t) (9) 

s 

where the sum is over the eight possible values of the three spins. The time dependence of the qk 
is given by 

qk{t) = -Mt) - ^[qk-lit) + Qk+lim (10) 
The solution to Eq. [10] is (see Appendix B ) 

qk{t) = ^fe(0)(e-(i-^)* + 2e-M)t] + b(fe-i)(0) + g,+i(0)] [e'^^"^)* - e~('~'^^'] (11) 

Evidently, at equilibrium, reached from any starting configuration, the average value of each 
spin is zero, as might be expected. 
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4.2 2-spin correlations 

The 2-spin correlation function is defined [Sj, in the 3-spin case, by 

rj,k{t) =< SjSk >= SjSkp{si, S2, S3; t) (12) 

s 

which guarantees that rjj = 1. The r's are, of course, symmetric in the two indices. The time 
dependences are given by 

7 

ri,i+i = -2rj,j+i + -[ri,i„i + ri_i,i+i] + 7 
7 

ri,j_i = -2rj,j_i + - [ri^i+i + rj_i,j+i] + 7 

7 

''j-i.i+i = -2ri_i,j+i + -[ri,i_i + rj,j+i] + 7 (13) 
The solutions are (see Appendix C) 

nAt) = r,,(0)e-^(^+^)* + ^i?(0)[e-(2-7)* _ e-(^(4+7)t] + _ ^-(2-7)*] (14) 

where R is the sum of the three independent correlation functions. The final term of Eq ll4lindicates 
that the 3 spins are partially correlated at all times. 



5 MAGNETIC FIELD 



Glauber observes that a dynamic Ising model with a magnetic field H will have the same equilibrium 
solution as the static Ising model if Wi in Eq [2] is replaced by 



where 



w'i{si) = Wi{si)[l - /3si] 



/3 = tanh(^) 



(15) 



(16) 



and /i is the magnetic moment associated with one of the spins. The Matrix M of the previous 
section is then replaced by 



-3(1-7) (1 + 7) 3(l-7)/3 (l + 7)/3 

3(1-7) -(1 + 7) -3(l-7)/3 -(l + 7)/3 

3(l-7)/3 (l + 7)/3 -3(1-7) (1 + 7) 

-3(l-7)/3 (3-7)/3 3(1-7) -(5 + 7) 



(17) 
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5.1 Zero Temperature: A Phase Transformation? 

The zero temperature ( 7 = /3 = 1) eigenvalues are readily found to be and 4, each occurring 



twice. The solutions may be written as 

xit) = {x(o) + 2t[x)o) -y(o)]}e-^* 
Y{t) = {y(o) + 2t[x)o)-y(o)]}e-^* 

U{t) = {[/(O) + X(0)(1 -e-^*) -2t[X(0) -y(0)]}e-^* 

V{t) = 1-Xit) (18) 

The equilibrium solution, which depends upon the initial conditions, has no counterpart in the 
original Ising model, where the concept of "initial" plays no role. The system ultimately relaxes to 
the configuration 

V{t) ~ 1 

U{t) ~ U{0) + X{0) 

= l-2V-{0) (19) 
where the last identity results from the fact that there are only three spins. 



A number of authors have treated the zero temperature limit of an Ising system as a phase 
boundary (see, generally, Privman in [4]). That this is a tempting analogy may be seen from the 
fact that in the zero temperature limit, the parameter |/3| is 1 when a magnetic field, no matter 
how small, is present, and if the field is absent. This parameter breaks the "parity" symmetry of 
the equilibrium solutions, resulting in the non-zero asymptotic values of Y{t) and U{t). 



5.2 Magnetic Field at Low, but Non-Zero Temperature 

The eigenvalues of the Master Equation matrix, Eq. [T71 must, in general, be determined from 
a cubic equation. I have investigated the solutions for the low, but non-zero, temperature case 
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where 7 = /3 = 0.90. I find the 3 non-zero eigenvalues to be approximately —4.81, —3.05, and 
—0.538. These may be compared with the values —6, —2.2. — 0.2, in the absence of a magnetic field 
{P = 0), and —4, —3.8, — .6[— 4, —2(1 + 7), —6(1 — 7)] in the presence of a very large magnetic field 
= 1). The magnetic field obviously breaks the parity symmetry, favoring +-spins over — spins 
(for positive IJ,H). 

The secular equation for the eigenvalues of the Master Equation matrix involve the parameter 
/3 to at least the second power, so that to linear order in the magnetic field the eigenvalues are 
unchanged from their zero-field values. To linear order in /3 the relaxation rates are just the zero- 
field rates. The solutions to the same order in /3, for a time dependent external field f3 = /Jqc^*'^* 
are: 

AU = -^([X(0)-34^]^^^^i±^(e-[2(2-7)+Ht_e-2(i-7)t) + 
2 + 7^^ ^ ^ 2(2-7)^ 2 + io; ^ ^ 

2(2 - 7 ^^_[2(2-7)+Ht _ e-6t) _ 3(1-7^) 2 + 7 ^^_^^^ _ ^_2(i-7)t^ 



Ay 



2(1 + 7 -iw)' ' (2(1 -7 -iw) 2 -7 

-J^/[X(0) - 3 (1-7) ] 47(l-7) ^_f2(2-7)+Mt _ g-2(l-7)t) _ 

2 + 7 2(2 — 7) 2 + iu! 

6(2-7) ^^-m2-y)+tw]t ^-6t^ 3(1-7)^ 2 + 7 ^^_^^^ ^_2(l-7H^^ ^20) 



2(1 + 7 -iw)' ' 2(1 -7) -iw] 2-7 

The equilibrium values are given by 

U(t) 



(4-7^)12(2-7)-H 

Y(f) 12/3o(l-7)^ iuii (2\) 

^ V^J ~ (4-72)[2(2-7)+ja;] ^'^'') 



Setting cj = makes it obvious that the magnetic field breaks the up-down symmetry of the no-field 
solutions. 

5.3 Average Values with a Weak Magnetic Field 

In place of Eq [10] the average value of spin k must satisfy, in the presence of a (weak) magnetic 
field H = Ho exp''^* 

Qk = -qk + 2^-1 + ^k+i) + /?[1 - -^{rk-i,k + rkMi)] (22) 

We suppose that the system is close to equilibrium, and approximate by using the equilibrium value 
2^ of the correlation functions rj.fc. The solution, after the transients have died off, is 

Qkit) ~ 3/3o(|^) . ^ (23) 

2 — 7 1 — 7 — ZCJ 

The average value of Glauber's stochastic magnetization", < M{t) >= fJ-J2kQk{t), is then 

<M(t)>=3i^^(|±I)-±^e-"' (24) 

kl 2 — 7 1 — J — loj 
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The time-delayed correlation function < M{0)M{t) > (see Appendix D)is 

< M(0)M(t) >= 3;u2(i±i) 6^(1-^)1*1 (25) 

^ 2 

The frequency-dependent susceptibility xi^)^ defined by < M(t) >= x{^)H is, then 

I ^ 2 + 7 1-7 

5.4 The Fluctuation-Dissipation Theorem 



The fluctuation-dissipation theorem, see e.g., |18) was probably motivated by Johnson's discovery 
of electrical "noise" voltage across otherwise isolated resistors |16j . 

Nyquist [17] subsequently related Johnson's voltage to "thermal agitation" of charges in the 
resistor, and predicted the magnitude of the voltage V in a frequency interval dv to depend upon 
the values of the ambient temperature T and the resistance R, according to 

V'^diy = ARkTdu (27) 

In the present context, following Glauber, we can let the magnetic field frequency uj simulate the 
frequency of the "thermal agitation" . Then the Fourier transform of the magnetization correlation 
function. Eg [25] is 

dt < M(0)M(t) > e-* = 3/i^(--^) 2 = -^Irnxiio) (28) 

Thus times the imaginary, or absorptive, part of the susceptibility, in this magnetic example, 
plays the role of the resistance in the case of Johnson noise. As Glauber points out, relations such 
as this may be useful in finding the effect of a weak field upon general functions of the spin variables 
in Ising-type models. 



6 Discussion and Conclusions 

The dynamic Ising model in one dimension would seem, at first blush, to be somewhat insipid — it 
can do nothing except relax to equilibrium. The paths to equilibrium, however, can involve all eight 
of the possible configurations, in the case of three spins. This raises the possibility of interesting 
patterns occurring in similar systems with larger numbers of spins. Even the three-spin model 
becomes more interesting, when used interactively, that is, by repeatedly shifting the equilibrium 
temperature. This is done by repeatedly changing the value of the 7 parameter, randomly or 
otherwise. 

The three-spin model with a weak magnetic field shows an amusing relationship to the infinite- 
spin model considered by Glauber. The susceptiblity in Glauber?s Eq. 94 (see also his Eq. 56) 
contains a term j^, which can be re-expressed as ^Jj^. The equivalent three-spin term in Eq[26l 
is just the small 7 (weak-coupling) approximation to Glauber's expression. 
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APPENDICES 

A TRANSITION PROBABILITIES 
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B AVERAGE VALUE SOLUTIONS 

Summing Eq. (jlOp over k (for arbitrary number of spins N), and defining Q{t) = J2kQk{'t) gives 
the equation 

Q(t) = -(l-7)Q(t) (29) 

with the solution 

Q{t) = Q(0)e-(i-^)* (30) 
Then, in the present instance of 3 spins we can rewrite Eq ([TO] ([14] Eq 32 et seq) 

q^{t) = -{l + l)qj,{t) + ^Q{t) (31) 

which has the solution 

q^it) = e-M)t[qUO) + Jq(0) /^^^dr] (32) 

^ JO 

which, in turn, gives Eqllll 

C AVERAGE CORRELATIION FUNCTION SOLUTION 

In keeping with the approach in Appendix B, we can define a quantity 

3 

R= ''i^k (33) 

i<k=l 
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The sum of the three equations ()13p then gives 

R- -{2--f)R + 2,-i (34) 

from which we deduce that 

R{t) = ii(0)e-(2-T)* + [1 - e-(2-7)i] (35) 

2-7 

but each of the equations (fT3|) mav be written as 

1 

ri,k = ^{R - n^k) + 1 (36) 
from which the solution given by Eg 1141 may readily be obtained by standard methods. 

D TIME DELAYED CORRELATION AND MAGNETIZATION 

The average two-spin time-delay correlation function 

< s,(0)sfe(t) >= r,j(0) exp-(i-^)* (37) 

is 



< Sj{0)sk{t) > 
1 
3 



< sj(0){sfc(0)[2exp-(i+i)* + exp-(i-^)*] 



sk-i{0) + Sk+m]} > [exp-(i-^)*-exp-(l + |)t]} 
1 



^,r,-fc(0)(2exp-(i+2)*+exp-(i"^)* - [r,- fc_i(0) + r,- fc+i(0)] 
(exp-(^+i)* -exp'(i-^)*) = njiO) exp-(i-^)* (38) 

: where the last equality reflects the fact that the rj^s are all equal at equilibrium. Summing over 
both indices (remembering that rjj = 1) and multiplying by fi'^ then leads directly to Eq. [25] 

E SIMULATION AND VIDEO DISPLAY 

There are eight possible configurations of the three two-valued spins. We number these as states 

= 1, .., 8 according two the scheme: 
N = 1(8) for configurations V+(y_) 
N = 1 + i for configurations x+j 
N = 3 + i for configurations 
The eight states are shown in eight .jpg files. 
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Transitions among the states, starting from some arbitrarily chosen initial state, with 7set 
in the range [-1,1], are encoded in an interactive python program which may be downloaded 
from http://www.hep.anl.gov/may/jlu/ThermoWith3spins/index.html. Just download and open 
the .tar file at that location and see the resulting README file for more information. 

The program simulates the relaxation to equilibrium of a 3-element Ising chain, and is repro- 
duced at the end of this Appendix, augmented with a graphics display. Programs, files and other 
supporting information can be found at the Web URL address cited above. The Web site also 
presents results of simulation runs as a graphical comparison of the simulation results with the a 
corresponding solution of the Master Equation. 

The graphical display of a run shows three arrows making transitions among the eight possible 
states. Each transition is one that is accomplished with a single-spin flip, as illustrated in the 
"cube" file, also available at the URL cited above. The transition rates, given in Appendix A, are 
assumed to be the rates of Poisson Processes [22], where the values in Appendix A are multiples of 
a basic rate "r" specified in The body of the program. 

The program starts with the python loop command "while k j= numberOflter" . As presently 
set up, this loop is first entered with the spins in the "ferromagnetic" state 1", displayed according 
the the variable "backl". Prior to the loop we define functions including eight functions denoted 
by "^n" (n=1..8), and a function "exponen(x)" . Each An function is called when the system is in 
the state corresponding to the value of "n" and returns three parameters on exit, namely, which 
of the three possible transitions will be made to a successor state, the rate at which the transition 
will be made, and the time spent in the current state prior to the transition. 

The Poisson flip rate for the successor state is given by the parameter "x" , specified in calling the 
exponen(x) function. That function randomly determines a waiting time z in the successor state 
in accordance with the exponential distribution of waiting times for the first event of a Poisson 
process. The probability, p(z), that the first event has occurred after elapsed time z 

p{z) = exp(— zx) (39) 

where p{z) is a random number in the range [0, 1] and z is the elapsed time variable (in the range 
[0,cx)]. Accordingly. 

z = —{In(randomnumberfrom[0, l]}/x (40) 

as has been emphasized by Gillespie) |21j . The variable z obviously represents the time spent in 
the successor state. 

During a single run of the simulation the fluctuations among the 8 possible states are large, 
making it difficult to "eyeball" an average time behavior to compare with the X(t), Y(t), ... func- 
tions given by the solution to the master equation. The conventional method is to use ensemble 
averaging |23j . by running the simulation using the same initial state but with different random 
number seeds. The states as a function of time are combined and sorted into small time bins for 
which the average values of X(t), Y(t),... are easily calculated. Typically results from 50 to 100 
runs are sufficient to obtain values in which the fluctuations are sufficiently small (i.e. standard 
deviation values for X(t) of 0.0806 and 0.0672 respectively) for reasonable comparisons to the mas- 
ter equation solutions. Both the master eequation relaxation times and asymptotic values are well 
reproduced by the simulation. Plots of these are available at the web site. 
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For a single interactive display run of the simulation which shows the 3-spin state as a function 
of time, we also show a plot of the X(t) and Y(t) function using a different averaging technique to 
give the user an indication of how the system relaxes to equilibrium. The algorithm accumulates, at 
each instant, the time spent in each of the states from time zero, thereby estimating the probabilities 
of each state as a function of time. This method reproduces the general behavior of the X(t), Y(t), 
and asymptotic values. 

A Python program to simulate the behavior of a 3 spin system using a Stochastic 

Monte Carlo Approach: 

# ! /usr/bin/python # this is a 
code to do a stochastic simulation of the dynamic time 
evolution of a state consisting of # 3 spin 1/2 objects 
based on a state transition diagram which specifies the 
available transitions # among the 8 possible states and the 
forward and backward transition rates. # Jack Uretsky and Ed 
May 2011 # this is glauberSspinSimulationS.py # this is 
revised to use exponential distribution in the A_ functions 
# to calculate the time of the next transition based on the 
rate local to the # A_ fimction # A graphical display (based 
on Tk functions) of the state and the X and Y functions as a 
function # of the iteration step. The display updates in the 
time determined by the simulation 

import sys import time import random import math 

debuglevel = 

def expon(x) : z = random . random ( ) exp Value = 
-math. log (z)/x return expValue 

if len(sys . argv) >= 1: print print "Usage: 
view3spinSimulation3 .py stateid gamma numberOf Iterations 
seed" print "For example: viewSspinSimulationS.py 1 0.25 325 
1379" time . sleep(2) # sleep for 2 seconds then continue 
with the simulation using default configuration #sys.exit() 
graphicsmode = "Tkinter" # viewSspinSimulationS.py [ Tkinter 
] stateid gamma numberOf Iterations runN seed # 
ViewSspinSimulationS.py stateid gamma 

numberOf Iterations seed # 1 

2 3 4 5 # decode the command 

arguement list to get initial parameters k=0 for arg in 
sys. argv: #print "Command line", arg, k k=k+l graphicsmode = 
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"Tkinter" from Tkinter import * 

initialstateid =1 # all spins Up if len(sys . argv) 

>= 2 : initialstateid = int (sys . argv [1] ) initialgamma = 0.25 
if len(sys . argv) >= 3 : initialgamma = float(sys.argv[2] ) 
numberOflter = 325 # fills the iteration length of the 

display 

if len (sys . argv) >= 4 : numberOflter = int (sys . argv [3] ) 
initialrunN = 1 
initialseed = 137 

if len(sys . argv) >= 5 : initialseed = int (sys . argv [4] ) 

if initialgamma <=-1.0 : 

print "Unphysical value for gamma. . .Exit" 
exitO 

if initialgamma >=+1.0 : 

print "Unphysical value for gamma. . .Exit" 
exitO 

print 

print "User requested graphicsmode=" .graphicsmode, "initialstateid=" .initialstateid 
print "gamma=" , initialgamma, "numberOf Iter=" .numberOflter, "initialseed=" .initialseed 
print 

from sys import exit 

import time 

import random 

from math import * 

global gamma 

gamma = initialgamma 

runN=initialrunN 

# setup structures to communicate with A_ functions 
global RET 

RET = [0,0,0] # will contain the triple state, rate, time given by the A_ functions 
global state 
global r 
global N 

# setup a dictionary to store state spin patterns for 3 spins 0=spin down l=spin up 
def bits (state) : 

s=' ' 
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#t={l:'lll', 2:'110', 3:'101', 4:'100', 5:'011\ 6:'010\ 7:'001', 8: '000'} #using apj 
t={l:'lll', 2:'011', 3:'101', 4:'110', 5:'100', 6:'010', 7:'001', 8: '000'} # jacks 21( 
s=t [state] 
return s 

# The A_ functions change the state at time t based on state diagram 
def A_l() : 

state =1 # current state 

r = 1 - gamma 

f = random. randrange (3) 

g = f+2 #destination 

t = expon( 3*r ) 

RET = [g,r,t] # g is new state after time t 

return RET 

def A_2(): 

state = 2 
global r2 
r2 = 1 . + gamma 
f = random. random 
if f > 2./ (3. + gamma) : 
g = 1 

elif f > l./(3. + gamma): 
g,r2 = 7,1 

else : 

g,r2 = 6,1 
t = expon(r2+2.0) 
RET = [g,r2,t] 
return RET 

def A_3() : 

state = 3 
global r3 
r3 = 1 . + gamma 
f = random. random () 
if f > 2./ (3. + gamma) : 
g = 1 

elif f > l./(3. + gamma): 
g,r3 = 7,1 

else : 

g,r3 = 5,1 
t = expon(r3+2.0) 
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RET = [g.rS.t] 
return RET 



def A_4(): 

state = 4 
global r4 
r4 = 1 . + gamma 
f = random. random 
if f > 2./ (3. + gamma) : 
g = 1 

elif f > l./(3. + gamma): 
g.r4 = 5,1 

else : 

g,r4 = 6,1 
t = expon(r4+2.0) 
RET = [g,r4,t] 
return RET 

def A_5() : 

state = 5 
global r5 
r5 = 1 . + gamma 
f = random. random 
if f > 2./ (3. + gamma) : 
g = 8 

elif f > l./(3. + gamma): 
g,r5 = 4,1 

else : 

g,r5 = 3,1 
t = expon(r5+2.0) 
RET = [g,r5,t] 
return RET 

def A_6(): 

state = 6 
global r6 
r6 = 1 . + gamma 
f = random. random 
if f > 2./ (3. + gamma) : 
g = 8 

elif f > l./(3. + gamma): 
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g,r6 = 4,1 

else : 

g,r6 = 2,1 
t = expoii(r6+2.0) 
RET = [g,r6,t] 
return RET 



state = 7 

global r7 

r7 = 1 . + gamma 

f = random. random 

if f > 2./ (3. + gamma) : 



g,r7 = 2,1 
t = expon(r7+2.0) 
RET = [g,r7,t] 

return RET 

def A_8() : 

state = 8 

global rS 

rS = 1 . - gamma 

f = random. randrange (3) 

g = f + 5 

t = expon(3*r8) 

RET = [g, r8, t] 

return RET 

#dictionaries 

Adestin = {1: A_l, 2 : A_2, 3 : A_3, 4: A_4,5 : A_5, 6 : A_6, 7: A_7, 8 : A_8} 
apics = {1: ' IMAGES/a_l . jpg' , 2: ' IMAGES/a_2 . jpg' , 3 : ' IMAGES/a_3 . jpg' , 4: ' IMAGES/a_4 . jpg' , 
5: 'IMAGES/a_5. jpg' , 6: ' IMAGES/a_6 . jpg' , 7 : ' IMAGES/a_7. jpg' , 8 : 'IMAGES/a_8. jpg'! 
state = initialstateid 
print "Beginning state=" , state 
iseed=initialseed 
random. seed(iseed) 

# open a data output file to store simulation data for this run 



def A_7() : 



g = 8 
elif f > l./(3 
g.r7 = 




else : 
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openf ile=open( 'glauberSspinSimulation.data' , 'w') 

# open a data output file to store configuration parameters for this rvm, used for plotting 
openf ileconf =open( 'glauberSspinSimulation. conf ' , ' w' ) 

backl = apics [state] # setup the initial state 

if graphicsmode == "Tkinter" : 

# Setup Tkinter graphics 

# create the canvas, size in pixels 

canvas = Canvas (width = 700, height = 750, bg = ' yellow O 

# pack the canvas into a frame/form 
canvas .pack (expand = YES, fill = BOTH) 

# create the X axis (ie iteration count) 
for i in range (50, 700) : 

canvas . create_text ( i , 625 , f ill= ' black ' , text= ' - ' ) 
asymX = (3*(l-gamina) )/(2*(2-gainma) ) 
j=625-100*asyiiiX 
for i in range (650, 700) : 

canvas. create_text(i,j ,f ill='black' ,text=' . ') 
canvas . create_text (600 , 625 , f ill= ' black ' , text= ' iteration' ) 
canvas. create_text(600, 525, fill='red' ,text='X() simulation') 
canvas . create_text (600 , 545 , f ill= ' blue ' , text= ' Y () simulation ' ) 

# create the Y axis at x=50 
for i in range (525, 725) : 

canvas. create_text (50, i,fill=' black' ,text=' . ') 

# create the Y axis tick marks 
for i in range (45, 55) : 

canvas. create_text(i, 525, fill='black' ,text=' . ') 

canvas . create_text ( i , 625 , f ill= ' black ' , text= ' - ' ) 

canvas. cr eat e_text(i, 725, fill=' black' ,text=' . ') 
canvas . create_text (30 , 525 ,fill=' black' ,text=' 1.0') 
canvas . create_text (30 , 625 , fill= 'black ' , text= ' 0.0') 
canvas . create_text (30 , 725 , f ill= ' black ' , text= ' -1 . ' ) 
c anvas . updat e ( ) 

k=0 

vPlusCnt=0 . 
vMinusCnt=0 . 
xPlusCnt=0. 
xMinusCnt=0 . 
pXsum=0 . 
pX=0 . 



18 



pY=0 . 

arbTime = 0. 

xAtTime =0 . 

xAtZero=0.0 

statePop= [] 

for i in range (8) : 

statePop . append (0) 



# write configuration data information 

openfile .write (' # iseed = '+str(iseed)+' init state = '+str (state)+' gamma = '+str(initialgami 
openfile .write ( '# numberOf Iterations = '+str (numberOf Iter) + '\n') 



# write out header 

openfile. write('# rowIndex'+' '+'state'+' arbTime ' '+'runN'+' '+'pX'+' '+'pY'+'\n') 

# write out starting state information 

openfile. write (str(k)+' '+str (state) +' '+str (arbTime) +' '+str(runN)+'\n') 



while k <= numberOflter : 
startTime = time. time () 
oldState=state 
k = k+1 

backl = apics [state] 

RET = Adestin [state] # this executes the appropriate A_ function to get a 

destin = RET[0] 
r = RET[1] 
state = destin 

if graphicsmode == "Tkinter" : 

# display the state using Tkinter graphics 
filename = 'IMAGES/a_'+str(state)+' .gif ' 
filename = 'IMAGES/a_'+str(oldState)+' .gif ' 

gif image = Photolmage (f ile = filename) 

canvas. create_image (20, 10, image = gif image, anchor = NW) 

j= 625-100*pX # calc a x,y position to plot pX value 

k2 = 2*k+50 

canvas . create_text (k2, j , f ill='red' , text='x' ) 

j= 625-100*pY # calc a x,y position to plot pY value 

canvas . create_text (k2 , j ,f ill='blue' ,text='+') 

c anvas . updat e ( ) 
# get time generated in the A_ functions 
rest = RET [2] 
arbTime=arbTime+rest 

StatePop [oldState-1] += rest # update a time weighted state 
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statelProb = float (statePop[0] )/arbTime 
# 

#t={l:'lll', 2:'011', 3:'101', 4:'110', 5:'100', 6:'010', 7:'001', 8:'000'} # jacks 21oct: 
#if (state==l) : vPlusCiit=vPlusCiit+l 
#if (state==8) : vMinusCnt=vMinusCnt+l 

#if ( state==2 or state==3 or state==4 ) : xPlusCnt=xPlusCnt+l 
#if ( state==5 or state==6 or state==7 ) : xMinusCnt=xMinusCnt+l 

# calculate various time weighted average for various X,Y,U,V functions at current time 
vPlusCnt = statePop[0] 
vMinusCnt= statePop [7] 

xPlusCnt = StatePop [l]+statePop [2] +statePop [3] 

xMinusCnt= statePop [4] +statePop [5] +statePop [6] 

pVplus = float (vPlusCnt)/arbTime 

pVminus = float (vMinusCnt)/arbTime 

pV= pVplus + pVminus 

pXplus = float (xPlusCnt)/arbTime 

pXminus = float (xMinusCnt)/arbTime 

pX = pXplus + pXminus 

pU = pVplus - pVminus 

pY = pXplus - pXminus 

#pXsum = pXsum + oldState*rest #? 

#pXave = pXsum/arbTime 

xAtTime=xAtZero*exp(-2* (2-gamma) *arbTime) +(3*(l-gamma)*(l .0 - exp(-2*(2-gamma)*arbTime)), 
#sys. stdout. write ('\r'+' at iteration='+str(k)+' state= '+bits(oldState)+' pX='+str (pX)+' 1 
print "\rAt iteration=y.d state=y.s pX=%6.4f time=y.6.3f" % (k,bits(oldState) .pX.arbTime) , 
sys . stdout . f lushO 

openfile. write (str(k) + ' '+str (state) + ' '+str(arbTime) + ' '+str(runN) + ' '+str(pX) + ' '+str(p'' 
time . sleep(rest) 

print "\n" 

vPlusCnt = StatePop [0] 
vMinusCnt= statePop [7] 

xPlusCnt = StatePop [l]+statePop [2] +statePop [3] 
xMinusCnt= statePop [4] +statePop [5] +statePop [6] 
pVplus = float (vPlusCnt)/arbTime 

pVminus = float (vMinusCnt)/arbTime 
pV= pVplus + pVminus 
pXplus = float (xPlusCnt)/arbTime 
pXminus = float (xMinusCnt)/arbTime 
pX = pXplus + pXminus 
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pU = pVplus - pVminus 
pY = pXplus - pXminus 

print "At end of iterations X=y.5.3f Y=y.5.3f V=y.5.3f U=y.5.3f" 7. (pX , pY , pV , pU) 
if ( debuglevel > ) : 
for i in range (8): 

print "y.5.3f " statePop [i-1] , 

print 

for i in range (8): 

print "y.5.3f " 1 (float (statePop [i-1] )/arbTime) , 

print 
print 

print "lastTime = y.8.3f" % arbTime 
print "gamma =", gamma 

print "Running at single spin flip rate (flips/sec) = %5.3f" % (numberOf Iter/arbTime) 
print "initialseed =" , initialseed 

print "Master Equation predicition at gamma = ", gamma, " and t — > infinity", 
print "X=y.5.3f" 1 ( (3*(l-gamma))/(2*(2-gamma)) ) 

openfileconf .write ("lastTime = "+str(arbTinie) + '\nO 
openf ileconf .write("gamma = "+str (gamma) \nO 
openfileconf .write ("initialseed = "+str (initialseed)+' \n' ) 
print "end of program" 



openf ile . close () 
openfileconf . close () 

print '\a', '\a', '\a' # ring terminal bell 
time . sleep(lO) 



print "That's All, Folks!" 
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